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0. Introduction /I * 

The objective of this report Is to present a new method of 
calculating viscous fluid flows having an average Reynolds number. 

The actual flows encountered In aeronautics have very large Reynolds 
numbers, but experience shows that the main flow structures are 
relatively Independent from the Reynolds number. 

In the short term, data extracted by a calculation using an 
average Reynolds number, such as the point of separation, the 
development and Importance of a separated flow, are of Interest 
to the aerodynamics specialist. In the long term, however, the 
numerical simulation of turbulent separated flows with large Reynolds 
numbers is searched for when refining an efficient and rellabj^ e Reynolds 
number, by adding to the Navler-Stokes equations either a standard 
turbulence model (TSEN-SEROU | l! ) or a homogeneous random turbulence 
(PERRIER-PIRONNEAU |2|). 

Furthermore, an effective numerical algorithm of 3-D unsteady 
dimensionless Navier-Stokes equations (1) is of Interest for 
Industry . 
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In (1) u designates the flow velocity, p the pressure, f the 
external forces, v the fluid viscosity (v » Re « Reynolds 
number), z is the velocity specification on the rigid or infinite 
boundary V of the Q volume occupied by the fluid, u^ the flow 
velocity at the intitial instant t » o. 

We propose the solution of (1) by the optimal control theory /2 
in terms of a system with distributed parameters (J-L LIONS | 3() 
and thereby to eliminate the two penalizing defects of the existing 
traditional methods: 

dl) Setting of parameters to ensure the convergence of the 
algorithm 

d2) Requirement of a very large computer. 

The decoupling method of the equations proposed (series of Stokes 
problems) is automatic (in the direction of dl) and requires a minimum 
band width (in the direction of d2). On the other hand, as it is more 
sophisticated, it shall be more difficult to program. 

The complexity of the aerodynamic geometries studied explains the 
preference given to a quantification of (1) into FINITE ELEMENTS . 

The large number of unknowns in Industrial anplications requires 
a rapid convergence (cr) for the treatment of non-linearities and 
schemes implicit in time (cr). 

I. QUANTIFICATION IN TIME AND TREATMENT OF NON LINEARITY /3 

To quantify the time derivative in (1), use is made of a scheme 
in two step s which is completely Implicit and unconditionally stable 
CCRANK-NICHOLSON) . If u”'^^(x) •* u(x, (n+1) At) deslpates the 
velocity at Instant (n+1) At, with At time step, u’^''’''’ is to be 
derived from u*^, u^“^ as the solution in (u,p) of (2) 
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with 


a • — tt f • ^3** _ 

‘ At At 


To solve the non-linear system (2), a least squares method in 
a Sobolev (fl), functional space, proposed in (GLOWINSKI- 
MANTEL-MERIAUX-PIRONNEAU j i| | ) , is given in (3) 


ma II (a- V A) S ♦ (u . V) u-f|j* 
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designates the dual space of 
(fl)** -{v c L* (m** / V V e (n)*^***. V Ir • o} 

K - 2 or 3 dimension of the region occupied by 


the fluid 
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where i ■ . 


to see that if (fi) contains a scalar product 

< u, V > - a (u, v) ♦ V (Vu, 7v) 

(5) 

• ) designates the f-calar oroduct in (G) 

(f, i) • /q H (6) 
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whereas (3) has the structure of an optimal 
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control problem (8 ) 
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in which J is the cost function or criterion and 


where u is a function of v via state equation (8) 
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TT acting as a pressure of which the variational formulation is (9) 


a jJjV, n d* ♦ V • j^£. n da - (v. V) V da 

v c I ¥ ^ C liS 


(9) 


(t) It may be noted that the state system is a modified Stokes 
problem. (GB4P|5j may be consulted to find the equivalence 
between the least squares — optimal control methods). 


II. SOLUTION OF AN OPTIMAL CONTROL PROBLEM BY THE CONJUGATE GRADIENT 
METHOD 


The solution of (7) and (8) is searched for the the POLAK- 
RIBIERE |6l version of the conjugate gradient algorithm. It exists 
is 3 steps: 


0. Initialization 


Select V® • (Stokes or idealized fluid) 

Calculate g° as the solution to the variational equation (10) 
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1. Descent (minimization to one variable) 


For m>o, assuming (v*”, as known, compute 


by 


X* - Arg sin J - > M*) I v», hP set 
X > o 


( 11 ) 


and where 




( 12 ) 


2 , Computation of a new descent direction 


Define g as the solution to the variational equation (1 

* ^ iP**» ^ d* ♦ V ^ Vn dK • « > ( 13 ) 

^♦1 e wg ; V M e 


Compute in (li4) / 6 _ 

>*♦' • a - g») dx ♦ V /pVj»*‘.V (J«Hl - i-j dx 

a /jj *** * '' I *** 

where ^♦1 . |j»»l ♦ Y^l > 

m ■ m+1 and refer to (11) 

Remarks: - (10) and (13) are Stokes equations 

-Minimization to one variable (11) Is achieved with the methofi 
of binary searches or with the Fibonacci (POLAK [6|) method. 
-The estimate requires detailing J* (v^**"^) 

Calculation of J’ and 

ij • < j*(v),6v • a /j^(u - V/ 4(u • v) dx 

♦ V /Qy(“ - v).t,*(Si - dx (16) 

♦ 0 (4u) ♦0(6^) 

where 6u is the solution of (17) via (8) 


5 


ro 



( 17 ) 
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Since <u-^) • Vq by using (16) (17) we have 
< J'(v),n > • a /g(u-v).n * v/jj7<u-v).7n ♦ /jj(u^).(v.V)r| 


n’ 


♦ /g (u-v) . 


¥ n fl «o 


( 18 ) 


To calculate from we must solve (8) 

with V ■ to obtain whereas by using (18) 

we h^ve > ■ 3 4j (u*** I ) .rd* ♦ v/qV(uP*i-v**') 

¥ntv, ^19) 

and is obtained from (13) (19) 
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In conclusion, each iteration of the altorithm (10).... (15) 
requires the solution of several Stokes problems: 

. The Stokes state eauation with v » to obtain 

. The gradient problem from u^'*’^, v”*'*’^ 

. The Stokes problems (- 3) for the evaluation of the cost 
function J during minimization to one variable (11). Consequently, 
an EFFICIENT amd RAPID Stoke? algorithm is a fundamental tool 
for solving the Navler-Stokes equations vlfi the least squares 
methods (3). 


III. METHOD OF DECOMPOSING THE STOKES ALGORITHP! 

(GLOWINSKI-PIRONNEAU l7| ) 

We limit ourselves here to the main concepts. For a more detailed 
development of the method, |8j may be consulted. The essence of the 


method consists In selecting the pressure on the boundary f in such a 
way that f »u ■ o Is satisfied In fl. 


Assuming the Stokes problem expressed In (20) 

<M - Au • ♦ f I V.tt • e I ttif. ■ u 


( 20 ) 


By taking the divergence of (20), we have 


a 


- 




( 21 ) 


If we know pU ■ X , we could compute the (iS,P) 
solution of (20) by solving (N+1) Dlrlchlet problems 

( 22 ) ( 23 ) 

- Ap • hi* in a ; pjj. - A ^22) 

a\ii • Auj - ll^ ♦ in n ; uilj. - i - l,M (23) 

Or the ^ solution ol (245 

M • ^*u In a j . o (24) 


By taking the Laplaclan ^ of (24) and by using (22) (24), 
we have 

-* « ■* -»♦ 

A A# • A (V*u) • V»(Au) • - Ap - y.f - aV.u - - oA# 


Or otherwise 

A A# ♦ oA# ■ o ; ♦|r • o 

Setting X so that || ]j. • o 

then * ♦ 

p i o St V.u • o 

Applying > — ^ U defined byreflnlng (21)-(24) 

we have A and b such that 


(26) 
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( 27 ) 




We must thus select X defined on the boundary as the solution 
to the linear problem (28) 




( 28 ) 


Remark: The solution of (28) is immediate, as A is syimnetrlcal and 12. 
highly elliptical . The functional support relating to these 
results is expanded in |9|- 

Construction of A 


The linear operator A is Implicitly defined by the sequence of 
Dlrlchlet problems (29) (30) (31) (32) 


® in * ^ 

(29) 

- -Vpj^ in « ^ “x*r * ° 

(30) 

- M) • ^*“a in ^ * * 

(31) 


(32) 


Construction of b 


b is explicitly defined by the sequence of Dlrlchlet problems 
(33) (3^0 (35) 


APo • in n ; Polj. • O 

oUq ■ auo • I - vpo i>» r. ; Uo'r • f 
- Mo • t.uo in n ; Solj. • o 


(33) 

(3^0 

(35) 
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(36) 


^^®(u,p) solution of (ao) is given by (37) 

( 37 ) 

with the trace of X of p|r as the only solution to the 
linear variational equation (38) 

< ax, u > - * ii. M > ® M defined in T 


IV, APPROXIMATION OP THE NAVIER-STOKES EQUATIONS BY THE FINITE /lO 
ELEMENTS METHOD 


If designates a polygonal approximation of the domain 
occupied by the fluid, is the set of triangles (Tj^) or TRIANGULATION , 
such that in a standard way : 

^ i TiA Tj - # If I # j 

h being tha largest side. 


The infinite functional spaces are substituted by the following 
finite spaces; 

K’ 1^- rolynomial with degree > ^ T stw/ 2 ' 

where is the trlangulatlcn obtained from ^ by dividing 
each triangl e into ^ sub~trlangles by drawing a line between the middle 
of the sides 


/;,i^ 
<^> (h.1, 
1^ .Uh • c« (r,) /g. 


Polynomial 

with 

degree 

2 ¥ T « C|, ) 

(1*0) 

polynomial 

with 

degree 

t V T iC|, ) 

(Hi) 

polynomial 

with 

degree 

1 V 1 S t!|, ) 

(Him 
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It 1 r stated that the approximation of the fluid velocity is 
linea r Cor quadratie). on the triangulation 

and that the approximation of the pressure is linear on the triangulation 



The problem is limited to discrete spaces (vj, P^, mJ) with 
the speed and pressure of Independent variables. 

If N2 designates the number of nodes of the triangulation /II 

N1 is the number of nodes of the trlangulatlon IT|, ant NIP the number 
of boundary nodes, whereas the number of unknowns of the discrete 
state system is NxN2+Nl+NlP with N - 2.3 dimension. 


In the sequence of equations of the Stokes algorithm, the linear 
solutions ■ d^^ are of three types: 

1) Discrete variational formulation of the state equation in 
velocity 


(^ 2 ) 




(4?) for \ . (no i - 1 , 52 , radix 

2) Discrete varational formulation of a pressure equation 




(^ 3 ) 


write (43) for #|| " (• using radix 

3) Discrete varational formulation of a pressure trace equation 

\iJ* * "I. • *s. ( w ) 


write 


(44) for • (w^) i • I, air, using radix > 1^1 


During an optimal control iteration, the number of dipcrete 
Stokes problems is approximately 5 


-one Stokes problem for the state equation 
-one Stokes problem for the calculation of the gradient 
of criterion J’ 

- 2^ 3 Stokes problems for the search for the optimal control 
Step w and for the direction h pre-deflned m^n J (w - Xh) 


Each Stokes problem requires (2xN 4) discrete scalar Dlrichlet 
problems as well as the solution of the small system (44), The empty 
Dlrichlet matrices corresponding to (- A) or (o Id -A) and the small 
complete matrice relatively to the pressure trace remain the same 
during the iterations and are thus factored out, once and for all 
(by a Cholevski direct type method) outside of the control 3oop. The 
construction of the small complete matrix (44) is a preliminary 
operation costing NlFx (N+2) Dlrichlet problems. 


At eacn time cycle nAt , several optimal control iterations (- 5) 
are required to ensure the convergence, the initial solution being the 
result of the preceding cycle (n-1) At. 

Numerical Simulations 


The optimal control method has boen t,.sted numerically on the 
I.E.M. 370/168. The How visualizations are the velocity field, the 
intensity of the rotational, the streamlines and the pressure fields 
in the case 2-D, the separated zones and the tube lines of the rotational 
for applications 3-D. 

"hp characterltties of each calculation (nodes, elements, Cholevski /I 
coefficient number, Reynolds number, calculation time C.P.U) are 
shown on the output plots of the code. The trlangulatlons are produced 
by the MODULE? techniques jl0|. The large number of s'^lutlons AX • B 
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to be found, with constant A, explains the preference given to a 
direct »CHOLEVSKI» or 'SKYLINE' method. The matrice band widths 
are minimized by the CUTHIL-MACKEE |ll| algorithms and A is factored 
out , once and for all, in a main core for simple test cases saved 
on secondary storage (disks) for complex industrial configurations. 

The Important role played by the Stokes solution (or of 
incompressible Idealized fluid) should be brought to light as a 
forecast for an internal flow estimate (or external) with a small 
Reynolds number. 

First of all, the validity of the code was tested on simple 
examples such as : 

1. The flow 2-D in a conduit with sudden enlargement to verify the 
convergence toward the steady state of the solution after several time 
cycles (Re = 100) Figures 1, 2, 3. Comparisons with the A. G. HUTTON 
code 1 12 I are presented. 

2. The unsteady 2-D flow around a Reynolds circle 200 with numerical 
perturbation after several time cycles (- 10) to produce behind 

the circle the Karmann path composed of alternating eddies Figure 

Finally, the flows on more complex geometries representing 
industrial configurations and requiring larger calculation times 
(> Ih CPU) were simulated. 

3. The separated flow 2-D around a profile of an air intake with a 
high incidence is analyzed at various time cycles. Figures 5 through 
13 show the origin, the development, the path and the disappearance 

of the eddies behind the profile (Re relative to the chord), on the 
Internal and upper external part of the air inlet (Re relating to the 
maximum deviation) and their effect on pressure. 

4. The separated flow 3-D behind a sphere and around a sweptback /lU 
wing with a high incidence. Figures 14 through 17 show the domain 

of return flow (u^ < o with U = (u^, , u^)) and several views of the 

tube lines of the rotational originating from the separated zone are 
12 presented • 
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NUMERICAL SIMULATION OF VISCOUS 
SEPARATED FLOW AROUND A IDEALIZED WING 


Reynolds 200 - Incidence 30* 
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